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The study of superconductivity with unconventional order is comphcated in condensed matter 
systems by their immense complexity. Optical lattices with their exceptional precision and control 
allow one to emulate superfluidity avoiding many of the complexities of condensed matter. A 
promising approach to realize unconventional superfluid order is to employ orbital degrees of freedom 
in higher Bloch bands. In recent work, indications were found that bosons condensed in the second 
band of an optical chequerboard lattice might exhibit px^ipy order. Here we present experiments, 
which provide strong evidence for the emergence of px Py order driven by the interaction in the 
local p-orbitals. We compare our observations with a multi-band Hubbard model and find excellent 
quantitative agreement. 

INTRODUCTION 

Understanding the role of unconventional order for superconductivity is a fundamental task in low temperature 
physics. Prominent examples in condensed matter are transition metal oxides [1 . Studying the order parameters in 
these systems is complicated by their vast complexity. A widely debated example is the chiral Px-\-'iPy order possibly 
formed in strontium ruthenates [2J| , which has recently attracted much interest because its topological nature may give 
rise to Majorana fermions \T . Optical lattices in their lowest band have proven to be a useful experimental arena to 
emulate superfluidity with exceptional precision and control [4^^. However, since under most general circumstances 
bosonic ground state wavefunctions are necessarily positive definite [6l[7] and hence topologically trivial, the realization 
of unconventional superfluid order in optical lattices with bosons is not straight forward. Possible approaches, presently 
receiving great attention, are based upon amending the scalar light-shift potentials of conventional optical lattices by 
static abelian and even non-abelian artificial gauge fields [8HTT] or upon use of dynamical lattice potentials p!2Hl4] . 
A method, more closely geared to electronic matter, is to employ atoms in metastable higher bands which provide 
orbital degrees of freedom [15 , 16 . Recently, we have shown that upon suitable control of band relaxation bosons can 
be condensed in the second band of a bipartite square lattice [17] . Momentum spectra were observed consistent with 
chiral Px^ipy superfluid order characterized by a spontaneously formed pattern of staggered local angular momenta, 
which breaks time-reversal symmetry. In the present work, by means of the following line of arguments we present 
clear evidence that Px^ipy is in fact formed: We show that for the lowest energy state of a bosonic Hubbard- model 
accounting for the second, third and fourth bands, the repulsive interaction in the p-orbitals stabilizes Px ^ipy order 
in analogy to Hund's second rule in multi-electronic atoms. We calculate a characteristic phase diagram with respect 
to a change of the interaction in the local p-orbitals and an adjustable distortion of the lattice, which tunes the energy 
minima of the second band. Experimental observations are presented, which show excellent quantitative agreement 
with the theoretical predictions of this phase diagram. We finally discuss excited state scenarios, which are compatible 
with the previously observed momentum spectra, but inconsistent with the experimental signatures reported in this 
work. 

LATTICE POTENTIAL 

We produce a two-dimensional optical potential comprising deep and shallow wells {A and B in Fig.[l](a)) arranged 
as the black and white fields of a chequerboard with an average well depth Vq and an adjustable relative potential 
energy offset AV pTHTO] . In the x?/-plane the optical potential is given by 

V{x, y) =-^\v (e^'" + ex e-'^-) + e'^ (e^^^ + e, e'^^y) \^ . (1) 
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FIG. 1: (a) The lattice possesses deep (A) and shallow (B) wells. A denotes the laser wavelength. A section of the potential 
along the dashed grey line is shown at the lower edge with the tunable relative potential energy offset AV indicated, (b) Contour 
plot of the second band within the 1st Brillouin zone for Vb = 6^rec, AV = 5.7 Erec- High and low energies are indicated 
by light grey and dark blue, respectively. The band exhibits two inequivalent local minima denoted as X± = |/iA:(l,±l) 
(k = A/27r) at the edge of the 1st Brillouin zone. 

Adjustment of /3 permits controlled tuning of AV = Vor]{l + ex){l + e^)cos(/3). A weak harmonic potential along 
the 2:-direction provides elongated tubular lattice sites. If 77 = e^^ = e^y = 1, the lattice potential possesses perfect 
C4 rotation symmetry. In our experiment, we are constrained to fixed parameter values r] = 1.03, and Cx = 0.93 
and hence C4 symmetry is weakly broken. In contrast to Ref.[17 , careful power and polarization management of the 
lattice beams permits controlled adjustment of arbitrary values of Cy around unity. The second Bloch band, shown 
in Fig. [l](b), provides two inequivalent local minima at the edge of the 1st Brillouin zone (denoted by X+ and 
which are energetically degenerate if the lattice displays C4 rotation symmetry. By adjusting the lattice distortion 
parameter Cy (see Appendix), we can continuously tune their energy difference AE = E{X_) — E{X_^) ^ 1 — Cy. In 
a tight-binding picture, the quantum states corresponding to X± may be approximated by Bloch states \^l^±) with 
real- valued Bloch functions V^±(r) composed of local p-orbitals {px, Py) in the deep wells and local s-orbitals in the 
shallow wells. Denoting their occupations per unit cell as Up and ris with no = + n^, the relative occupations 
i/p = Up/no and = rig /no only depend on the spatial shape of V^±(r) (but not on the local chemical potential) and 
hence can be tuned via adjustment of AV (see Appendix). 

MULTI-BAND HUBBARD MODEL 

In order to predict the nature of the expected quantum phases for different values of Up and AE^ we employ a bosonic 
multi-band Hubbard-Hamiltonian H (see Appendix) accounting for the three tight-binding bands associated with the 
three local orbitals Px^ Py and s. We include nearest-neighbour (NN) and next-nearest-neighbour (NNN) tunneling 
processes, on-site interactions for the p- and 5-orbitals, and a term accounting for a potential energy difference between 
p- and 5-orbitals, required for modeling the tuning of AV. Furthermore, an extra term is included, which introduces a 
quadrupolar anisotropy of the tunneling between p-orbitals with respect to the {cx + /i e2y)-directions (with Cx^y shown 
in Fig. [l](a) and /i G { — 1, 1}). This term with the real amplitude Jy^^ permits to model the tuning of AE according 
to the relation AE ^ 8upJ\\^a^ which is obtained by evaluating the single-particle spectrum of H in momentum space 
at X± . Comparing the lowest two tight-binding bands of H with the second and third bands of a full two-dimensional 
band calculation for the experimentally realized lattice potential lets us determine the tunneling parameters. The 
largest tunneling amplitude Jsp arises for NN-tunneling between s- and p-orbitals (see Appendix). 

For AE = 0, any linear combination |6>,0) = sin(^) + cos{0) e^^\ip-) minimizes the single-particle energy. For 
repulsive collisions, minimization of the interaction energy in the p-orbitals requires maximal angular momentum 
[151 [16] ^iid hence (j) = ±7r/2 and = ii/A because the local superposition states Px ± ipy are maximally delocalized 
such that the atoms can optimally avoid each other. For AE ^ 0, depending on the sign of A£^, either |?/^+) or |?/^_) 
minimize the single-particle energy. When AE is increased beyond some critical value, the gain of single particle 
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FIG. 2: The phase diagram in the centre shows the mixing angle versus rip 



(in units of Jsp)- Three phases arise 



(regions (I), (II), and (III)) separated by 2nd order transitions (white dashed hues). On the upper edge the corresponding orbital 
and phase ordering is illustrated. Orbital currents and plaquette currents are highlighted by circular and straight arrows. The 
colors and numbers indicate the local phases of the different local orbit als. Below the phase diagram momentum spectra are 
shown for different values of J\\ „. 



energy exceeds the cost of interaction energy required for eliminating angular momentum, and hence also the total 
energy is minimized by one of the states \^|^±). A mean field analysis for the general case AE ^ shows that the ground 
state of R can be in fact approximated as |^,±7r/2) = sin(6>) \^^) ± icos{6) |?/^_) with the mixing angle 6 plotted 
in the phase diagram in the centre of Fig. [2] versus Jy^^i and UpUp {Up is the on-site collision energy per particle in 
the Px- and p^-orbitals). Note that sin^(^) and cos^(^) quantify the relative populations of the condensation points 
X+ and X_, respectively. This phase diagram comprises three regions, separated by second-order phase boundaries 
highlighted by white dashed lines. In regions (I) and (III) one finds ^ = and = 7t/2 respectively (hence, only one 
of the condensation points X± in Fig. [l](b) is occupied), while in region (II) the simple relation 

cos(2.) = -i^ = -^ (2) 
^ ^ ripUp 2nlUp ^ ^ 

holds [20] . The mixing of the two condensates is governed by a competition between the gain of single-particle energy 
per unit cell no AE introduced by the lattice distortion and the interaction energy per unit cell Up gained by 
maximizing angular momentum in the p-orbitals. The pictograms on the upper edge of the phase diagram illustrate 
the orbital and local phase ordering predicted in the three regions. In regions (I) and (III) the order parameters 
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are real and their local phases indicated by the colors and numbers are arranged in order to maximize tunneling, 
while interaction energy does not play a role. In region (II) an inherently complex- valued order arises with orbital 
currents and plaquette currents highlighted by circular and straight arrows, respectively. Interaction energy sets the 
relative local phases between p^- and p^-orbitals at the same site to be 7r/2, while the phase relations between orbitals 
at neighboring sites are arranged to maximize tunneling. The unit cell of the order parameter comprises four unit 
cells of the lattice potential and time-reversal is equivalent to a shift by one unit cell of the lattice potential. Hence, 
time-reversal symmetry is broken. Our theoretical considerations are consistent with a numerical analysis based upon 
the Gross-Pitaevskii equation fST and a renormalization group analysis [22^ . 

EXPERIMENT 

Our experimental procedure begins with a Bose-Einstein condensate of rubidium atoms (^''Rb) loaded into the 
lowest band. By means described in Ref. [TT the atoms are excited to the second band. Their temperature remains 
sufficiently low, such that after typically 10 ms they condense in the energy minima of the band. In our experiments 
we wait for 80 ms to ensure complete thermalisation. Below the phase diagram in Fig. [2j we show momentum 
spectra observed in the different areas (I), (II) and (III). In regions (I) and (III) the observed Bragg-resonances 
unambiguously prove the realization of the predicted standing-wave order in configuration space sketched in the 
corresponding pictograms. In region (II), where both points X+ and X_ are populated, the momentum spectra 
appear as superpositions of those in regions (I) and (III). This clearly confirms that the underlying quantum state 
is composed of two states \^|J±) corresponding to condensates at X± with Px±y order. However, a more precise 
determination of its nature requires additional information. In the following, we show that the state prepared in our 
experiment closely follows the phase diagram in Fig. [2] derived for the ground state |6>, ±7r/2) of which we refer to 
as scenario A. 

In Fig. [3] we compare the observed populations of the condensation points X± as a function of Up and the lattice 
distortion with the predictions of scenario A. The experimental procedure yielding the data points (filled red disks) in 
Figs. 3 (a)-(f) is as follows: For different settings of the lattice distortion parameter e^, AV is varied adiabatically and 
the normalized mean occupation difference between both condensation points (z^dif) = ((^+) ~ (^-))/((^+) + (^-)) 
is recorded and plotted versus Up. Here, n± are the populations observed in X± in a single measurement and (...) 
denotes the average over multiple experimental realizations. For the small interaction energies realized in our system 
- with Up/Jsp on the order of a few percent - and temperatures well below the condensation temperature, our mean 
field zero temperature analysis appears well adapted to model the observations if the finite size of the lattice with 
spatially varying local values of Up and Up is accounted for. In fact, upon applying a local density approximation, 
the observations are remarkably well reproduced by the theoretical predictions for scenario A. Using according to 
Eq. I^jwith the local values of AE^ Up and Up^ we calculate the corresponding local value of z^dif,th = sin^(^) — cos^(6>) 
at each plaquette in the lattice and subsequently apply an average over all plaquettes to obtain I^dif,th7 which is 
plotted as the solid green lines in Figs. 3 (a)-(f). The grey areas refiect the largest uncertainty in this calculation given 
by the calibration of the lattice distortion parameter. A discussion of the technical details is given in the Appendix. 
According to the observations in Fig.[3j n+ and n_ approach each other if the distortion is reduced or if Vp is increased, 
in excellent quantitative agreement with the predictions for region (II) of the phase diagram in Fig. [2j where the value 
of J|| and hence AE required to produce imbalanced condensate fractions grows with increasing Vp. This striking 
behavior is hence identified as a signature of Py :^ipy-OTdei described by scenario A. The physical mechanism may be 
regarded as an analogue of Hund's second rule for repulsively interacting bosons: maximal local angular momentum 
is favorable in order to minimize interaction energy. 

Aside from the ground state scenario A, the distribution of the Bragg-resonances in region (II) is also compatible 
with three possible excited-state scenarios: the coexistence of spatially separated real phases \ip±) (scenario B), a 
real coherent superposition (1 T l)'7r/2) = sin(^) |?/^+) ± cos(^) |?/^_) with Px+y Px-y order (scenario C), and 
an incoherent mixture described by the density operator po = sin^(^) |?/^+)(?/^+| + cos^(^) |?/^_)(?/^_ | (scenario D). 
The compelling quantitative agreement of our observations with the predictions for the ground state scenario A is 
underlined by the fact that the excited state scenarios B, C, and D lead to predictions in explicit contrast to the 
observations. The phase separation scenario B exhibits spatially separate Px-\-y and Pa^_^-orbitals corresponding to the 
two condensates \^l^±). Hence, interaction energy does not favor equal populations of the condensation points and thus 
^dif,th cannot acquire any dependence upon Up^ in sharp contrast to the observations. The absence of phase separation 
in our experiment is also supported by the fact that in our Bragg spectra we never observe different coherence areas for 
the Bragg peaks corresponding to the two condensation points. This would be expected, if the measured condensate 
populations n+ and n_ would be attributed to spatially separated domains of different size. The irrelevance 
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FIG. 3: The normalized mean occupation difference (z^dif) between the two condensation points X± is plotted versus the 
relative occupation of the p-orbitals Vp. The error bars show the statistical errors for eight measurements. The adjusted 
distortion of the lattice potential is 1 - = 0.019, 0.014, 0.009, 0.004, 0.001, -0.001 from (a) to (f). The solid green lines show 
the corresponding theoretical predictions X^dif,th derived by means ofEq.[2] The grey areas represent the uncertainty of the 
lattice distortion Ae^ = ±2.5 • 10~^. The blue dashed lines show the predictions for scenario D. 



of scenario B in our experiments is not surprising: Coexisting spatially separated domains of real phases |?/^+) and 
|?/^_) as in scenario B are energetically more costly than either of the pure phases IV^i). A change of the domain sizes 
merely requires a local redistribution of particles between the Px+y and p^^-^-orbitals. Particle transport over many 
lattice sites is not necessary. Hence, equilibration should occur within a few tunneling times leading to the formation 
of either |?/^+) or IV^-). The real superposition state of scenario C exhibits spatially separate Px and p^-orbitals. 
For equal populations of the condensation points every second local s-orbital remains unoccupied due to destructive 
interference, which is energetically strongly unfavorable. In fact, minimization of the energy of|^,(l^l)7r/2) with 
respect to 9 yields either of the values 6> = 0or^ = 7r/2 only depending on the sign of and hence i^dif,th can only 
take the values ±1 (see Sec. X of Appendix). Finally, although the incoherent mixture pe in scenario D provides 
orthogonal Pa;±2y"0^t)itals at the same lattice site, the interaction energy gained by equally distributing the atoms 
among these orbitals is smaller than for scenario A because of the indeterminate phase relation between the two 
states IV^i). We have minimized the energy of pQ with respect to in order to determine I^dif,th for this scenario. 
The result (see Appendix, Eq.A3) is plotted as the dashed blue lines in Fig. [sj which obviously disagrees with the 
observations. In fact, significantly larger values of Vp would be required to equilibrate the condensate fractions as 
compared to scenario A. Scenario D also appears implausible because the two incoherently superimposed condensates 
can exchange particles through binary collisions in the shared local s-orbitals of the shallow wells (see pictograms of 
regions (I) and (HI)). This should rapidly degrade the coherence in each condensate, which is not observed. 

Further insight is gained by analyzing the fluctuations of v± = n±/{n-^ + n_). Histograms recorded for ^ 
show fluctuations of — well described by Gaussians, centered at U-^ = V- (see Fig. [4] (a)). This observation 
appears incompatible with the phase separation scenario B, for which, similarly as in Ref.[14 , a double peak structure 
should be expected instead, since equal populations n± of the two involved condensation points are not preferred. In 
Fig.[4](b) the standard deviations Ai^sum = (((^+ + ^-)^) - (^+ + ^-)^)^^^ and Ai^dif = (((i^+-i^-)^)-(i^+-i^-)^)^/^are 
plotted versus v^. The observations show that — = {Au^^^ — Au'^^^) /A significantly deviates from zero, 

i.e., the fluctuations of u± are strongly correlated. While the comparably small fluctuations of the total condensate 
fraction Az^sum show basically no dependence on Az^dif is sizable and notably decreases as Up is increased, showing 
that the interaction in the p-orbitals tend to lock the populations condensed in X+ and X_ . This observation again 
clearly rules out scenarios B and C, which do not prefer equal values of n± as Up is increased. It clearly reflects the 
physics of scenario A captured in the phase diagram in Fig. [2] As Up is decreased, the critical point is approached 
along the vertical Jy = 0-line where all three phases adjoin. The larger Up the higher the interaction energy to be 
paid for deviations of U-^ — V- from zero yielding suppression of fluctuations Av^a. 
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FIG. 4: (a) Histograms of zy+ — v- for Vp — 0.06, 0.10, and 0.32 each showing more than 150 identical reahzations. The red 
soUd traces show Gaussian fits, (b) The standard deviations Az^sum and Az^dif are plotted versus Vp (open and filled blue disks: 
Az^dif, black squares: Azygum; the open disks correspond to the histograms shown in (a); the thin grey lines connecting the data 
points are for guiding the eye). 



CONCLUSION 



Optical lattices with ^ ipy order open exciting perspectives for future research. Matter wave interference 
techniques [23 could be used to further study the mutual coherence of the condensates at X+ and X_. Imaging 
of the atoms with single-site resolution as demonstrated in Refs. [24l [25] might allow one to directly observe the 
local angular momentum of the wave function and hence to explore the spontaneous symmetry breaking process. 
Proceeding to deeper potential wells, one may access the strongly correlated regime where a rich phase diagram 
of Mott insulators with distinct orbital ordering is expected [26-28 . One may also explore topologically protected 
features in higher bands [29 and, if fermions are used, simulate forms of topological matter [30 resembling those 
discussed in the context of electronic systems. 



APPENDIX 



A I. Lattice potential. Using an interferometric lattice set-up [18], we produce a two-dimensional optical potential 
comprising deep and shallow wells {A and B in Fig. [l] (a) of the main text) arranged as the black and white fields of 
a chequerboard with an average well depth Vq and an adjustable relative potential energy offset [IZ^. In the x7/-plane 
the optical potential is given by 

U(x, y) =-^\v {e'^'' + e-^^") + e'^ (e^^^ + Cy e'^^y) p . (3) 

Adjustment of /3 with a precision exceeding 7r/300 permits controlled tuning of AV = Vq r]{l -h ^x){^ + ^y) cos(/3). A 
weak harmonic potential (with 40 Hz vibrational frequency) along the z-direction provides elongated tubular lattice 
sites. If 7^ = = e^y = 1, the lattice potential possesses perfect C4 rotation symmetry. In our experiment, due to 
unavoidable imperfections of the lattice set-up, we are constrained to fixed parameter values 77 = 1.03, and Cx = 0.93 
and hence C4 symmetry is weakly broken. In contrast to Ref.[T7l, the optical set-up permits controlled adjustment 
of arbitrary values of Cy within an interval including Cy = 1. This is accomplished as follows: the optical standing 
wave along the 7/-axis is obtained by a retro-reflected laser beam. The linear polarization of the incoming beam can 
be rotated with a retardation plate. After retro-reflection the polarization is rotated to precisely match with the 
2;-direction, which exclusively contributes to the lattice potential. 

A II. Band structure: tight-binding picture and full-band calculation. One may understand the 
single-particle bands within a simplified tight-binding (TB) picture. As sketched in Fig. |5](a), each local vibrational 
orbital of the two types of wells gives rise to a Bloch band. We operate in the regime where the 2nd, 3rd and 4th 
bands are significantly closer to each other than their separation from the 1st band, which correspond to the local 
5-ground states in the deep wells. We may write AU = Vsp + AUgp with a potential energy offset Vsp depending on 
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FIG. 5: (a) Schematic of first four bands plotted in configuration space for AVsp < (upper detail) and AVsp > (lower 
detail). The black and red dashed lines represent degenerate local px- and p^-orbitals yielding two closely spaced bands. In 
(b) and (c) the 2nd (red), 3rd (green) and 4th (blue) Bloch bands are plotted along a roundtrip in the 1st Brillouin zone (grey 
square in (d)) connecting the points F, X+, M, X_, F. (b) and (c) show the idealized C4-symmetric case (?7 = = = 1) and 
the experimentally implemented case (77 = 1.03, = 0.93, = 1), respectively, (e) TB bands based upon the Hamiltonian 
in Eq. Q fitted to match the bands in (b). (f) Tunneling amplitudes used in the Hubbard model in Eq. Q. The numbers u 
plotted inside the local s- and p-orbitals denote their local phases e**^^"^-*^/^. 



the adjusted value of Vb, and AVsp denoting the potential energy difference between the local 5-states in the shallow 
wells and the local Px and Py orbitals in the deep wells (cf. Fig. [s] (a)). For Vq = 6£^rec used in our experiment, 
Vsp = 5.7£^rec- If AVsp < O5 the second band arises primarily from the local 5-states in the shallow wells (upper 
graph in (a)), while for AVsp > it corresponds to a superposition of the degenerate Px- and p^-orbitals in the deep 
wells (lower graph in a). A two-dimensional band calculation for the potential in Eq. (|3| yields the true 2nd, 3rd and 
4th Bloch bands. In Fig.[5](b) and (c), these bands are plotted along a roundtrip in the 1st Brillouin zone connecting 
the points F,X+,M, X_,F (indicated in (d)). The average well depth is Vq = 6£^rec and AV = 5.7£^rec, which 
corresponds to AVsp ^ in the TB picture in (a). In (b) the idealized C4 symmetric case is shown, compared in (c) 
to the experimentally implemented case with Cy = 1. The band degeneracies arising at the F- and M-points in (b) are 
lifted in (c) due to the lattice distortion. In either case the 2nd band possesses two inequivalent band minima at the 
X_- and X+-points with energies E{X±)^ which are degenerate even in the case (c) despite the broken C4-symmetry. 
Only if Cy is tuned away from unity, this degeneracy is lifted. Experimentally, tuning of AE = E{X-) — £^(X+) is 
accomplished by tuning of Cy with AE{ey) ^ l — Cy quantified by a band calculation for the potential in Eq. ([3|. In the 
vicinity of the X_- and X+-points, where the condensed atoms reside, the bands in (b) and (c) are approximately equal. 

A III. Bosonic multi-band Hubbard model. We employ the multi-band Hubbard-Hamiltonian 

H = -Jss ^R^R+f^d^ - Jsp ^ (MPa,i?^^+/ie. + ^.c) 

Yl iJ\\^^J\\,a)pi,RPa,R+f,d^ -J± Yl ^ ( P^RPx^R+f^d^) (4) 

ReA,(j,fi,u ReA,fi,u 

+ Y - ^p,r) + ^ ^^l,R ~ ^1,r/^) ^ Y ^ rLs,R{^s,R - 1) , 

ReA ReA ReB 

with the summation indices jJi^i' ^ { — 1,1}, cr G {x,?/}, dy = Cx ^ i^ey^ and Cx^y as shown in Fig. [l] (a) of the 
main text. This Hamiltonian accounts for all possible tunneling processes between nearest- (NN) and next-nearest 
neighbors (NNN), as is indicated in Fig. |5](f): NN-tunneling between 5-orbitals and p-orbitals (J^p), NNN-tunneling 
between 5-orbitals (Jss), NNN-tunneling between p^^-orbitals or p^-orbitals (Jy), NNN-tunneling between px- and 
p^y-orbitals ( J±), a term accounting for a potential energy difference between p- and s-orbitals (AVsp) and the on-site 
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interactions for p-orbitals (Up) and 5-orbitals (^/s), respectively, with hs^R = Sj^sr^ hp^R = p\. j^Px^R Py RPy,R 

and Lp^R = i{pl RPy,R — p\ rPx.r)- Furthermore, an extra term scaling with Jy is included, which introduces a 
quadrupolar anisotropy of the tunneling between p-orbitals with respect to the d^y- direct ions. This term permits to 
model the tuning of the energy difference between the minima of the second band. Diagonalizing the kinetic part 
of H in momentum space at X± yields the relation /S.E ^ SvpJ\\^a with the relative occupation Up of the p-orbitals 
given as i.p « i[l + AVsp/{32J^p + AV;pi/2]. 

A IV. Determination of hopping parameters and phase diagram. The Hubbard Hamiltonian of Eq. Q 
yields three TB bands associated with the three local orbitals involved. Setting Jy = we match the lowest two TB 
bands with the full bands in Fig. [5](b). This lets us determine the tunneling parameters as Jss ^ 0, Jgp = 0.12£^rec5 
J|| = 0.07 Jgp, J± = 0.15 Jgp, and AVsp = 0.3 Jgp. The resulting TB bands are plotted in Fig. [s] (e) showing good 
agreement with the full bands. The smaller bandwidth of the 3rd TB band as compared to the corresponding 4th 
full band indicates the proximity of higher bands (p-orbitals in the shallow wells) in the true lattice potential, not 
accounted for in the TB-description. A non-zero value of Jy ^ introduces an imbalance for NNN-tunneling between 
p-orbitals along the x-\-y and x — ^/-directions, which acts to lift the degeneracy of the two band minima of the lowest 
TB band. The phase diagram in Eq. (1) of the main text is derived as follows: the Hamiltonian in Eq. Q including 
J||^a is rewritten, using a composition of the lattice via plaquettes with four sites. We diagonalize the kinetic energy 
in momentum space and assume that the system is condensed at the two inequivalent energy minima of the lowest 
TB band. Applying a mean-field approximation, we calculate the interaction and minimize the total energy with 
respect to the relative phase between the two possible condensation states and with respect to their relative weights. 

A V. Momentum spectra and band populations. Momentum spectra are obtained by rapidly (< 1 /is) 
extinguishing the lattice potential, permitting a free expansion of the atomic sample during 30 ms, and subsequently 
recording an absorption image. Band populations are measured as follows: the population of the n-th band is 
transferred into the n-th Brillouin zone by adiabatically terminating the lattice potential in 400 /is, followed by a 
ballistic expansion of 30 ms. An absorption image of the atomic density distribution is recorded and the populations 
in the different Brillouin zones are counted. 

A VI. Tuning of Up^ effect of band relaxation. The required tuning of the relative population of the 
p-orbitals Up is accomplished via adjustment of AV. Since this tuning is essential, we have studied it in some detail. 
The value of Up corresponding to some /3 and Cy (and hence AV and AE) is obtained by integrating |V^±(r)p over 
those fields of the chequerboard lattice comprising the deep wells with the Bloch functions V^±(r) at the condensation 
points X± derived from a band calculation employing the potential of Eq. ([3|. Vp scales monotonously with AV (as 
shown by the solid black trace of Fig. [g] (a)) and it is practically independent of the lattice distortion parameter Cy 
and hence of AE. This behavior is consistent with the approximate analytic expression derived from the Hubbard 
model in Sec. A HI. 

Changing A V also has a significant impact on the time-scale of band relaxation, as is shown by the data points in 
Fig. [6](a). The band lifetime is found to be maximal (~ 230 ms) if most atoms reside in the local 5-orbitals {up ~ 0). 
This is expected, since for these atoms there is no local state with lower energy available, which could give rise to 
relaxation. The initial preparation of the condensate is carried out in this configuration. Following a holding time 
of 80 ms to reach complete equilibrium, a subsequent adiabatic increase of Vp during 3 ms is directly visible in the 
momentum spectra: the higher order Bragg-peaks become more populated due to the increased contribution of the 
local p-orbitals, which due to their nodal structure comprise higher momenta. This is shown in Fig. [g] (b) and (c) 
with the corresponding calculation in Fig.[6](d) yielding good agreement. Increasing values of Up are accompanied by 
faster band decay. Since a single tunneling process between adjacent wells is sufficient for adjusting Up^ adiabaticity 
only requires tuning times of a few tunneling times of about 1 ms. 

The lifetime of the atoms in the second band (Fig. |6] (a) red disks) is measured as follows: after forming the 
condensate at AV = 3£^rec7 is tuned in 3 ms (sufficiently slow to permit tunneling) to the desired value. After a 
variable duration the population of the second band is determined (see V.). The wings of the decaying populations 
are fitted by exponentials with the 1/e-times plotted as the red disks in Fig. [g] (a). The lifetimes of the condensed 
fraction (blue squares in Fig. [6](a)) are obtained by an analog procedure, however counting the number of atoms in 
the lowest order Bragg peaks of a momentum spectrum. 

A VII. Calibrating the lattice distortion Cy. Changes of AE are experimentally implemented by ad- 
justment of the lattice distortion parameter Cy using polarization optics. Calibration of AE{ey) proceeds as follows: 
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FIG. 6: (a) The data points show the hfetime of the atoms (red disks) and that of the condensed fraction (blue squares) in 
the second band versus AV for Vb = 6.0 ^rec- The soUd hne (derived from a band calculation) shows the relative p-orbital 
population Up versus AV. The vertical dashed-dotted line indicates AV — Vsp = 5.7 Erec, where z/p = z/^ = 1/2. (b) Two 
observations of momentum spectra recorded for small (3.0 ^rec) and large AV (6.5 ^rec)- The red dashed rectangle identifies 
selected Bragg-maxima analyzed in more detail in (c) (observations) and (d) (theory) for varying values of AV^. 



the linear polarization of the incoming beam in the ^/-branch of the lattice potential is rotated such that (n+) ~ (^-), 
which corresponds to AE{ey = 1) = 0. Arbitrary values of AE are adjusted by rotating the polarization away from 
this position by precisely quantified amounts and determining the corresponding values of AE via a band calculation 
for the potential in Eq. ([3|. The estimated error in the determination of Cy is Acy ~ 2.5 • 10~^. 



A VIII. Determination of UpUp^no AE and I^dif,th' We account for the isotropic harmonic potential with 
approximately 40 Hz trap frequency superimposed upon the two dimensional lattice of Eq. (|3| with r] = Cx = Cy = 1 
and tubular lattice sites extending along the z-direction. The total number of condensed atoms, determined by fitting 
Gaussians to all visible Bragg peaks of a momentum spectrum and counting the atoms, is ^ 1.5 • 10^ before band 
relaxation sets in. Assuming a Thomas-Fermi density distribution along the weakly confined 2;-direction, we calculate 
the local number of atoms per unit cell tiq^r at the Bravais lattice site and hence the local number of particles in 
the s-orbitals and p-orbitals, Us^r = i^s'^o,R and Up^R = VpUo^R^ respectively. In the center of the lattice uq^r ~ 50. 
In order to determine the local interaction energies per particle in the p-orbitals and s-orbitals at site Up^R and 
Us^R^ the local wells at the A- and S-sites are approximated to 6th order and the corresponding wave- functions of 
the local p-orbitals and 5-orbitals are calculated. With this input, we determine the local energies ^ Up^R and 
riQ^R AE. By applying Eq. 2 of the main text, the corresponding local value of and hence i^dif,th = sin^(^) — cos^(^) 
is obtained. Averaging over all positions in the lattice yields n'^Up^noAE^ used in Fig. ^ (a) and I^dif,th plotted in 
Fig. [3] of the main text. 

A IX. Measurement of (vdif) ^ind comparison to theory i^dif,th- In order to obtain Fig. [3] (a)-(f) of 
the main text, in Eq. ^ Cy is adjusted to the six different values 1 - G {0.019,0.014,0.009,0.004,0.001,-0.001} 
and P is varied for each setting. Momentum spectra are recorded and the number of atoms n± in the two zero 
order Bragg peaks corresponding to X± are counted. A data point in Fig. [s] (a)-(f) represents an average over 8 
measurements. We thus obtain (i^dif) = ((^+) — (^-))/((^+) + (^-)) ^^r different combinations of the parameters 
P and Cy. From /3 one obtains AV (see Sec. A I) and thus Up (see Sec. A VI). The value AE is determined from a 
band calculation for the potential in Eq. ([3| (see Sec. A VII). For each data point, no and Up = VpUo are derived at 
the time of the measurement, thus accounting for band relaxation loss. With this the energies Up^riQ AE and the 
occupation difference i^dif,th are obtained (see Sec. A VIII). The data points in each of the graphs (a)-(f) in Fig. [s] 
of the main text correspond to a path in the plane spanned by the energies Up and no AE. These pathes are 
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FIG. 7: (a) Map of the values of rip Up and no for the data points in the graphs (a)-(f) in Fig. pi of the main text. The 
kinks in the paths arise due to coUisional losses setting in for large values of Vp. (b) For the case of the red dashed path (b) 
corresponding to Fig.[3](b), the local value of z^dif,th(i?) is plotted across the lattice. The lattice is parametrized by the Bravais 
vectors R = fid^ + v d-. 



identified here in Fig. [7](a). For the case of the red dashed path corresponding to Fig. |3](b), the local value of i^dif,th 
is plotted across the lattice in Fig. [t] (b). This confirms that, particularly for the data points in the center of the 
path, finite size effects are in fact substantial and have to be accounted for. 

A X. Mean field predictions of I^dif,th for real superposition state and incoherent mixture. We 

have also minimized the total energy for the real coherent superposition sin(^) |?/^+) ±cos(6>) |?/^_) (scenario C) and for 
the incoherent mixture of spatially superimposed condensates sin^(^) \ +cos^(^) |?/^_)(?/^_| (scenario D) in or- 

der to calculate the local mixing angle 6> as a function of and and hence to determine i^dif,th = sin^ (O) — cos^ (0) . 
For scenario C we obtain the simple result that i^dif,th is constrained to the values ±1 depending on the sign of A^ 
with no dependence upon v^. For scenario D, values of i^dif,th deviating from ±1 require Vp > 2/3 and satisfy 

3A^ 

[TJ^ -4(1 - lypHnoUp 

An average over all plaquettes of the lattice as decribed in Sec. A VIII leads to the corresponding value of i^dif,th, which 
yields the (blue) dashed traces in Fig. [s] of the main text. 
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